%% "Econometrics of the Hodrick-Prescott Filter," Robert de Jong and Neslihan Sakarya (2015)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This file allows you to compare Theorem 1 representation of the HP filter with the HP filter code written by Ivailo Izvorski, the Department of Economics,  Yale University.
% The link for the code: https://dge.repec.org/codes/izvorski/hpfilter.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% After you run the program, the following results are displayed in the command window:
%
% (1)A (Tx1) vector in which the sum of the weights are calculated for each t \in {1,2,...,T}. Each entry of the vector should equal to 1 (sum_{s=1}^{T} w_{Tts}=1)
 
% (2)A (Tx4) table that allows you to compare 4 series (in order)
%    (i) The original series; y (a random walk process is generated) [labeled in the table as "Data"]
%   (ii) The trend of y that is calculated by using Izvorski's code [labeled in the table as "HPfilter"]
%  (iii) The trend of y that is calculated by using the weights of the HP filter given in the Theorem 1 [labeled in the table as "HPfilter(Thm1)"]
%   (iv) The difference between the two trend series [labeled in the table as "HPfilter-HPfilter(Thm1)"] 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Explanation of the codes below:
% In this program, we first write the functions "eta, phi, f_T_lambda, and g_T_lambda"
% to calculate the weights of the HP filter as represented in Theorem 1 of the paper. 
% Then, we give a code for the HP filter written by Ivailo Izvorski. Izvorski's code allows us to see the accuracy of the HP filter representation in Theorem 1.
% Lastly, we generate a random walk series {y_t} and detrend it in two
% different ways; (1) using Izvorski's code and (2)using weights of the HP
% filter given in Theorem 1  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
clear all;
close all;
clc;
T=180; %% total number of observations
lm=1600 ; %% the smoothing parameter, lambda

%% Defining the Building Block Equations: delta(dlt), eta, xi, phi
dlt=0; 
eta=0;

for j=1:T;
    
    dlt=dlt+((sin(pi*(j-1)/(2*T))^2*cos(pi*(j-1)/(2*T)))^2*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1)); 
    eta=eta+((sin(pi*(j-1)/(2*T))^4*cos(pi*(j-1)/(2*T)))*cos(pi*(j-1)*(T-(1/2))/T)*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1));
end

dlt=(1/T)*dlt; %delta
eta=(1/T)*eta; %eta
ksi=32*lm*(1-64*lm*dlt)*(1-64*lm*dlt+(32*lm)^2*(dlt^2-eta^2))^(-1)+(32*lm)^2*(1-64*lm*dlt+(32*lm)^2*(dlt^2-eta^2))^(-1)*dlt; %xi
phi=(32*lm)^2*eta/(1-64*lm*dlt+(32*lm)^2*(dlt^2-eta^2)); %phi

%% Defining the function f_T_lambda, it is denoted as f

a=zeros(T,T); % the matrix for cos(pi*(j-1)*m/T) // used in f_T_lambda function 
df=zeros(T,1); % the vector for 1/(1+16*lm*sin(pi*(j-1)/(2*T))^4) // used in both f_T_lambda and g_T_lambda

for j=2:T;
    for m=1:2*T;
        a(m,j)=cos(pi*(j-1)*m/T);
        df(j)=1/(1+16*lm*sin(pi*(j-1)/(2*T))^4);
    end
end
b=a*df; %% sum j=2 to T cos(pi*(j-1)*m/T)*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1)

f=zeros(2*T,1); % the maximum value of m should exceed sample size T, since m might take larger values than T i.e. when m=t+s-1 and take t=T and s=T-1 then m=2T-2 > T

for m=1:2*T; 
   f(m)=(2*T)^(-1)+cos(pi*m)*(2*T)^(-1)*(1+16*lm)^(-1)+T^(-1)*b(m); %% f function (Note that cos(pi*m)=(-1)^m)
end

%% In this part, f_T_lambda (0) is defined since there is no zero entry of the vector f, it is denoted as f0 and all the corresponding functions will be denoted similarly, 
 
% Note that cos(pi*(j-1)*m/T)=1 when m=0

b0=0;
for j=2:T;
    b0=b0+(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1); % sum j=2 to T (1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1)
end

f0=(2*T)^(-1)+(2*T)^(-1)*(1+16*lm)^(-1)+T^(-1)*b0; % Note that cos(pi*m)=1 when m=0


%% Defining the function g_T_lambda, it is denoted as g

r=zeros(T,T);
dg=zeros(T,1);

for m=1:T;
    for j=1:T;
        
        r(m,j)=sqrt(2)*cos(pi*(j-1)*(m-(1/2))/T); % the matrix for sqrt(2)*cos(pi*(j-1)*(m-(1/2))/T), used in g function 
        dg(j)=sin(pi*(j-1)/(2*T))^2*cos(pi*(j-1)/(2*T))*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1); %(q1(j)*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1)), used in g function
    end
end

g=T^(-1)*r*dg; %% r*dg is the sum j=1 to T (sqrt(2)*cos(pi*(j-1)*(m-(1/2))/T)*q1(j)*(1+16*lm*sin(pi*(j-1)/(2*T))^4)^(-1))

%% Defining TxT weight matrix: w_{T,t,s} is the (t,s) entry of weight matrix with the number of observations T
 %Defining Indicator functions I(t+s-1=T)=I1, I(t+s-1<T)=I2, and I(t+s-1>T)=I
for t=1:T;
    for s=1:T;
        
       if t+s-1==T;
           I1(t+s-1)=1; %% indicator function I(t+s-1=T)
       else
           I1(t+s-1)=0;
           
           
           if t+s-1 < T;
               I2(t+s-1)=1; %% indicator function I(t+s-1<T)
           else
               I2(t+s-1)=0;
               
               
               if t+s-1 > T;
                   I(t+s-1)=1; %% indicator function I(t+s-1>T)
               else 
                   I(t+s-1)=0;
               
               end
           end
           
       end
    end
end

wts=zeros(T,T); % the weight matrix
f0ts=0;

    for t=1:T;
        for s=1:T;
            if abs(t-s)==0;
                f0ts=f0; %% f0ts is a function that takes value f0 when t=s, f(abs(t-s)) otherwise
            else 
                f0ts=f(abs(t-s)); %% f(abs(t-s)) is taken instead of f(t-s) since f(t-s)=f(abs(t-s)) by symmetry of cosine function
            
            end
            %the formula of the weight of point s for the point t
            wts(t,s)=f0ts+f(T)*I1(t+s-1)+f(abs(t+s-1))*I2(t+s-1)+f(2*T-t-s+1)*I(t+s-1)+ksi*g(t)*g(s)+phi*g(T-t+1)*g(s)+phi*g(t)*g(T-s+1)+ksi*g(T-t+1)*g(T-s+1);
        end
    end
    
    


%% Showing the equivalence of the HP filter formula (using Izvorski's code) and Theorem 1 representation of the HP filter 

  %Generating a random walk series y and the error term e
  
  y=zeros(T,1);
  e=randn(T,1);
  y(1)=randn(1,1);

    for j=1:T-1;
    
    y(j+1)=y(j)+e(j+1);
end
 
 
%The HP filter code written by Ivailo Izvorski.

  if size(y,1)<size(y,2)
   y=y';
  end
  
t=size(y,1);
z1=6*lm+1;
z2=-4*lm;
z3=lm;
d1=[z3,z2,z1];
d1=ones(t,1)*d1;
m1=diag(d1(:,3))+diag(d1(1:t-1,2),1)+diag(d1(1:t-1,2),-1);
m1=m1+diag(d1(1:t-2,1),2)+diag(d1(1:t-2,1),-2);

m1(1,1)=1+lm;       m1(1,2)=-2*lm;
m1(2,1)=-2*lm;      m1(2,2)=5*lm+1;
m1(t-1,t-1)=5*lm+1; m1(t-1,t)=-2*lm;
m1(t,t-1)=-2*lm;    m1(t,t)=1+lm;

trend1=inv(m1)*y; %trend1 is the trend of y series calculated by original HP filter formula

% The HP filter representation in Theorem 1
 
 trend2=wts*y; %multiplying weight matrix with the original series y gives the trended y series

%% The results that are displayed on the command window

% the sum of the weights of a particular point t should be equal to 1      
  sumw=wts*ones(T,1);
  printmat(sumw, 'sum of the weights')
  
% A table that compares 4 series; y (a random walk process), the trend of y calculated by Izvorski's code (trend1), the trend of y calculated by the representation in the Theorem 1 (trend2),
% and the difference between two approaches (trend1-trend2)

R=[y, trend1, trend2, trend1-trend2];

disp('   Data                HPfilter           HPfilter(Thm1)         HPfilter-HPfilter(Thm1)')
disp(R)
